function [ tree, rates ] = buildEMTree( copyNums, headers )
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
    maxCopyNum = max(max(copyNums)) + 1;
    rates = randInitRates(maxCopyNum);
    thresh = .0001;
    converged = 0;
    iter = 0;
    while (converged == 0)
        % Get estimated copy number profiles for internal nodes
        dists = getDistsWithTransitionProbabilities(copyNums, rates);
        tree = seqneighjoin(dists, 'equivar', headers);
        allCopyNums = getAllCopyNumsRound(tree, copyNums);
        
        % Get new rates given internal nodes
        ratesNew = getTransitionProbs(tree, allCopyNums);
        
        dif = sum(sum(abs(rates - ratesNew)));
        fprintf('at iter %d, diff is %d\n', iter, dif);
        diag(ratesNew)
        if (dif < thresh)
            converged = 1;
        end
        rates = ratesNew;
        iter = iter + 1;
    end
end

function rates = randInitRates(maxNum)
    rates = rand(maxNum);
    for i = 1:maxNum
        rowTotal = sum(rates(i, :));
        rates(i, :) = rates(i, :) / rowTotal;
        if (abs(sum(rates(i, :)) - 1) > .000001)
            fprintf('prob error! found %d\n', sum(rates(i, :)));
        end
    end
end

